Transition to chaos and modal structure of magnetized Taylor–Couette flow

Taylor–Couette flow (TCF) is often used as a simplified model for complex rotating flows in the interior of stars and accretion discs. The flow dynamics in these objects is influenced by magnetic fields. For example, quasi-Keplerian flows in Taylor–Couette geometry become unstable to a travelling or standing wave in an external magnetic field if the fluid is conducting; there is an instability even when the flow is hydrodynamically stable. This magnetorotational instability leads to the development of chaotic states and, eventually, turbulence, when the cylinder rotation is sufficiently fast. The transition to turbulence in this flow can be complex, with the coexistence of parameter regions with spatio-temporal chaos and regions with quasi-periodic behaviour, involving one or two additional modulating frequencies. Although the unstable modes of a periodic flow can be identified with Floquet analysis, here we adopt a more flexible equation-free data-driven approach. We analyse the data from the transition to chaos in the magnetized TCF and identify the flow structures related to the modulating frequencies with dynamic mode decomposition; this method is based on approximating nonlinear dynamics with a linear infinite-dimensional Koopman operator. With the use of these structures, one can construct a nonlinear reduced model for the transition. This article is part of the theme issue ‘Taylor–Couette and related flows on the centennial of Taylor’s seminal Philosophical Transactions paper (part 1)’.

AG, 0000-0003-2831-184X Taylor-Couette flow (TCF) is often used as a simplified model for complex rotating flows in the interior of stars and accretion discs. The flow dynamics in these objects is influenced by magnetic fields. For example, quasi-Keplerian flows in Taylor-Couette geometry become unstable to a travelling or standing wave in an external magnetic field if the fluid is conducting; there is an instability even when the flow is hydrodynamically stable. This magnetorotational instability leads to the development of chaotic states and, eventually, turbulence, when the cylinder rotation is sufficiently fast. The transition to turbulence in this flow can be complex, with the coexistence of parameter regions with spatio-temporal chaos and regions with quasi-periodic behaviour, involving one or two additional modulating frequencies. Although the unstable modes of a periodic flow can be identified with Floquet analysis, here we adopt a more flexible equation-free data-driven approach. We analyse the data from the transition to chaos in the magnetized TCF and identify the flow structures related to the modulating frequencies with dynamic mode decomposition; this method is based on approximating nonlinear dynamics with a linear infinite-dimensional Koopman operator. With the use of these structures, one can construct a nonlinear reduced model for the transition.
This article is part of the theme issue ' Taylor

Introduction
Stability and transition to turbulence in fluid flows have been of interest for scientists since the beginning of the twentieth century. In 1923, Taylor explored a linearly unstable flow between two concentric rotating cylinders both theoretically and experimentally in his influential work [1]. Its starting point was a combination of the Rayleigh stability criterion for inviscid rotating fluids, for instability, (1.1) together with the viscosity measurements by Couette and Mallock, indicating an instability in an analogous set-up. Here, angular velocities Ω i , Ω o and radii r i , r o correspond to the inner and outer cylinders. Taylor's stability diagrams for axisymmetric disturbances are in excellent agreement with experiments though the final pages of his work are devoted to the instability of the axisymmetric vortices themselves. He observed that 'a large increase [in the speed] caused the symmetric motion to break down into some kind of turbulent motion· · · ', and that 'each vortex was pulsating so that its cross-section varied periodically'. Since then, Taylor-Couette flow (TCF) has also become a model for turbulence generation in rapidly rotating astrophysical flows [2,3], where turbulence is important for angular momentum transport, magnetic field generation and mixing of chemical species. In those studies, the velocity of the cylinders is set to approximate the desired astrophysical rotation, for example, Keplerian profile Ω ∼ r −1.5 of an accretion disc. Quasi-Keplerian flows, with their angular momentum increasing with radius, are hydrodynamically linearly stable to infinitesimal perturbations according to the Rayleigh criterion (1.1). Though transition to turbulence through finite perturbations at large rotation speeds cannot be ruled out, TCF is remarkably stable in quasi-Keplerian regimes, up to Reynolds numbers of Re ∼ 10 6 in experiments [4]. Thus, other physical mechanisms of instability in quasi-Keplerian flows are frequently considered. One of them, magnetorotational instability (MRI), arises in differentially rotating flows threaded by large-scale magnetic fields, frequent in astrophysical objects. TCF was used as a model for experimental studies of MRI [5].
The stability of TCF in the presence of magnetic fields was first studied by Velikhov [6] and Chandrasekhar [7] in the 1950s. Considering axisymmetric perturbations and axial magnetic field, Velikhov concluded that magnetized TCF is unstable if i.e. when the angular velocity, and not angular momentum, decreases with radius. In an ideally conductive fluid, a radially displaced fluid parcel drags away the magnetic field line, 'glued' into the flow and retains its previous angular velocity. In the new location, the fluid parcel experiences three forces: magnetic tension, the centrifugal force and the equilibrium pressure gradient. If the velocity decreases outwards, the centrifugal force of the fluid element is larger than the pressure gradient, and in sufficiently weak fields this leads to instability. If the magnetic field is too strong, magnetic tension stabilizes the flow. Similar arguments can be invoked for azimuthal magnetic field; here the flow stability depends on the radial shape of the field [6,8]. Hollerbach et al. [9] showed numerically that non-axisymmetric disturbances with azimuthal wavenumber m = 1 are the most unstable in this case. Most of the existing MRI studies focused on asymptotic behaviour of instability and properties of fully developed turbulence. Transition from MRI to turbulence in TCF was not investigated systematically. Guseva et al. [10] found that MRI arises in a supercritical Hopf bifurcation, and then the flow undergoes a subcritical Hopf bifurcation to chaos when the ratio of viscosity ν to magnetic diffusivity η of the fluid is low, i.e. magnetic Prandtl number Pm = ν/η ∼ 10 −6 . In plasma-like fluids with large Pm ∼ 1, Guseva et al. [11] reported a more complex scenario of transition, with the flow passing a succession of oscillatory states as the strength of magnetic field varies. The present work aims to analyse how these oscillatory states appear and evolve using a data-driven approach. We employ the method of dynamic mode decomposition (DMD) for the identification of coherent structures in physical systems. DMD was developed by Schmid [12] as an alternative to costly iterative methods of global stability analysis, with direct applications to fluid flows. It can be interpreted as the generalization of global stability analysis for both numerical and experimental data, and results in a set of 'dynamic' spatial modes and corresponding eigenvalues. We will refer to them as DMD modes and DMD eigenvalues, respectively. For nonlinear systems, DMD represents linear tangent approximation of the system's dominant dynamics. The theoretical significance of this linear approximation is closely related to the idea that the dynamics of a nonlinear system of finite dimensions can be represented by a linear infinite-dimensional Koopman operator [13]. This operator propagates flow observables in time, and its eigenvalues and eigenvectors fully define the dynamics of the system; DMD can be viewed as the numerical approximation of this operator. Compared to other decomposition methods like principal orthogonal decomposition (POD), DMD is superior in identifying flow frequencies, and therefore is more appropriate for the analysis of the above-mentioned oscillatory states.
This paper is structured as follows: first, we introduce our numerical set-up for TCF, and the qualitative description of the transition to turbulence. After that, we describe the DMD method in more detail. We present the results of DMD analysis and identify the dynamical components related to the transition. Finally, we discuss the results and give an outlook on the possible future work.

Description of the flow
The equations describing the motion of an incompressible conducting fluid in the presence of magnetic fields are the Navier-Stokes and induction equations and Note the feedback of magnetic Lorentz force (∇ × B) × B on the flow; this force is an essential component of MRI. Thus, the flow is intrinsically nonlinear in velocity field u and magnetic field B. The laminar solution to (2.1) in hydrodynamic TCF is The radius ratio of the cylinders was set to r i /r o = 0.5, and the rotation rate to Ω o /Ω i = 0.26, approximating a quasi-Keplerian profile with Ω → r −2 , but still fulfilling criterion (1.2 implies that the dissipation of magnetic and velocity fluctuations takes place on the same scale, and leaves only two free parameters in the flow, Re and Ha. We solve (2.1)-(2.3) using direct numerical simulations (DNS) in Taylor-Couette geometry. The laminar velocity profile (2.4), and the azimuthal magnetic field B φ = B 0 (1/r) are imposed as forcing terms. The code has spectral discretization in the axial and azimuthal directions z and ϕ, and the radial coordinate r is discretized using finite differences. The nonlinear terms are evaluated in the physical space and are de-aliased using the 3/2 rule; more details of the numerical method can be found in [10]. The axial wavenumbers are set to k = αk , k = 0, 1, . . ., with α = 0.5, equivalent to setting the length of the cylinders to L z = 2π/α = 4π . The spatial resolution is N r , N z , N ϕ = (40, 64, 16), where N z and N ϕ in the number of Fourier modes in respective directions.  Each run is started from a small non-axisymmetric perturbation of one flow mode with k = 2 (αk = 1), m = 1. This initial condition simplifies the flow dynamics, constraining it to the subspace with only even wavenumbers, k = 0, 2, 4, . . .. Without this constraint, the odd wavenumbers also become active, which results in a considerably more complex transition scenario left for the future work.

(a) Magnetorotational instability and transition to chaos in direct numerical simulations
The linear stability analysis of the flow, performed by linearizing (2.1), (2.2) and with the numerical method from [9], shows the parameters Re, Ha where MRI is active (figure 1a). The instability arises when Re > Re cr ≈ 100; above this threshold, the magnetic field should be neither weak nor too strong for the flow to be unstable. As the Reynolds number increases, the magnetic field strength required to trigger instability also increases; however, the instability range becomes wider overall. We focus on the case with Re = 250, as in [11]. The instability growth rates ω r along this line are represented by the dashed line in figure 1b.  In DNS, the diagnostic quantity for the onset of instability is the friction torque on the cylinders G. It is related to transverse current of azimuthal motion in radial direction J ω [14] through where the angular momentum can be transported through the tension of magnetic field lines via the Maxwell stress component B r B ϕ [3]. Here · · · A denotes a spatial average along a cylindrical surface A. J ω is constant across the radius r, so that G i = G o = G on average for statistically steady flows. In the absence of MRI, the laminar flow torque G lam is constant in time and can be calculated analytically from (2.4). As the instability develops, the friction on the cylinders increases. This increase is directly related to the dissipation enhancement in the flow [15], as more energy is required to maintain the rotation. Figure 1b shows a correlation between the increase in G and the instability growth rate ω r . Both G and ω r reach their maximum at about Ha = 120, however, G is not monotonic, with a local minimum developing where ω r is the largest. Finally, the instability ceases to exist at about Ha = 230 as the magnetic tension becomes too strong and stability is restored. Now we focus our attention on the temporal behaviour of G. The instability appears first as a standing wave at Ha ≈ 50 (figure 2a), which corresponds to a time-independent friction and energy state once initial transients saturate (figure 1c). With increase in Ha, the flow rapidly becomes chaotic, with an abrupt transition to chaos at low Ha and only a narrow interval of doubly periodic in G solutions. At Ha = 100, the velocity and magnetic fields exhibit chaotic features (figure 1g), but retain some spatial structure (figure 2d,e). The vertical component of velocity and magnetic field is periodically dominated by a large-scale structure (figure 2f ). As Ha is increased further, the chaotic behaviour begins to regularize, and a slow modulation becomes discernible at Ha = 120 (figure 1f ). Soon, the torque timeseries is nearly non-chaotic and again shows doubly periodic behaviour (figure 1e), with a rapid oscillation of ω/Ω i ≈ 0.09, and the slower modulation of ω/Ω i ≈ 0.02. Neither of these frequencies corresponds to the frequency of the MRI mode which rotates much faster azimuthally, at ω/Ω i ∈ (0.  also reduces with further increase in Ha, until it ceases to exist at Ha = 149, and the steady state is again a standing wave ( figure 1c). The strong asymmetry in transition to chaos (abrupt on the left, gradual on the right) is possibly related to subcriticality of the left stability border of the MRI [3], although further work is necessary to confirm this. Table 1 gives an overview of the Ha intervals of the transition from chaos to regular behaviour. In the following, we will focus on this parameter region, decreasing Ha from about Ha = 150 to Ha = 100, so that the flow complexity increases.  figure 3h. Overall, the flow transition to chaos through the breakdown of a torus falls into the Ruelle-Takens scenario of transition to turbulence [16,17]. Nevertheless, here even chaotic flow states retain some regularity, and therefore could be potentially described with a few relevant dynamical components. In the next section, we will use the data-driven method of DMD to approach this problem.

Dynamic mode decomposition
Consider the system (2.1) and (2.2) in a general form, (u, B)  (e)  where f is a nonlinear operator. We seek the best linear approximation to this nonlinear system in the form of In general, the eigenvalues of (3.2) are complex, i.e. ω = ω r + iω i . In the simulations, the information about the flow is available in the form of magnetic and flow field snapshots q, sampled every t in time, so it is more practical to seek a discrete-time system (i) Collect snapshots q k of the system at timesteps k = 1, 2, . . . , K.
Step 4 introduces the truncation parameter r, typically defined by some criterion in the spectrum of singular values σ in the diagonal of matrix Σ in 3. By increasing this parameter, more dynamical information about the system may be kept. However, as singular values in decomposition 3 decrease, the singular vectors associated with these small singular values are increasingly linearly dependent and including them in 4 would make the subsequent decomposition ill-conditioned [12]. A common choice is a 99% cut-off of the SVD spectrum σ . In our case, this criterion leads to a very large intractable dynamical basis for the flow. Since the singular values spectrum σ 2 represent the energy content of the POD modes U, we re-define the cut-off parameter r so that r σ 2 r / σ 2 = 99%, with the modes retaining 99% of the energy of the respective quantity. Different components of magnetic and flow field have different cut-off values, depending on their spatial complexity. In general, u r , B r have the largest modal basis, and u ϕ , B ϕ , influenced heavily by their mean fields, remain low-dimensional. For the latter, we set the cut-off parameter at r = 3, including the mode with ω i = 0, corresponding to non-oscillating motion, and two main complex-conjugate frequencies. The optimal amplitudes of each mode, indicating its relative importance for the flow, were calculated as a best-fit of the data onto the DMD model [19].

Results
In the following, we discuss our DMD results, corresponding to the dynamical regimes in table 1. As the flow variables are related through (2.3) and nonlinear terms in (2.1) and (2.2), they have similar frequency content, so we first focus on the axial velocity u z . Figure 4 presents its DMD spectra in the form of discrete eigenvalues (3.3). Figure 5 depicts the spatial structure of some of the modes of u z and B z , and figure 6 compares DMD results to our linear stability analysis and DNS. We begin by discussing first the common features of the DMD spectra at different Ha, and then focus on the transition between the flow regimes.

(a) Magnetorotational instability and non-oscillatory modes
In the DMD spectrum, the main oscillating component of the signal is represented by two modes with complex-conjugate frequencies, denoted by 1 and 1 * in figure 4. This is the dominant mode for the explored values of Ha, and we denote its continuous-time analogue as ω 1 . Its dominant wavenumber is k = 3, with six pairs of rolls along the domain length L z = 4π (figure 5a); it is nonaxisymmetric with azimuthal wavenumber m = 1, like the standing wave flow pattern in figure 2a.
In figure 6a, we compare the dominant DMD frequencies of u ϕ and B ϕ , (ω 1 )/Ω i ∈ (0.3, 0.4), with the frequencies of the dominant MRI wave from the linear analysis. This frequency represents the azimuthal rotation of the MRI wave; both the linear and DMD frequencies decrease with Ha, with the latter slightly higher up until Ha ≈ 120. Then, as the emerging chaos becomes more pronounced, the dominant DMD frequencies become smaller than the linear ones, especially for the modes of magnetic field. This frequency adjustment is expected, as the nonlinear saturation of the instability modulates its initial growth, and the chaotic flow includes several dynamical components with comparable frequency content (figure 4d). Despite the emerging chaotic motion, the agreement between the frequencies indicates that the MRI-unstable modes remain active. The rest of the modes in figure 4 are clustered around ω 1 and its higher harmonics, ±nω 1 , n = 2, 3  MRI modes in the nonlinear terms of (2.1), (2.2), and have a finer spatial structure of (k, m) = (6, 2). They do not represent independent dynamics and are unlikely to play a role in the transition to chaos here, since it involves frequencies slower than the MRI ones. The non-oscillatory modes with (ω 0 ) = 0, corresponding to a purely real discrete eigenvalue λ = 1, are denoted by 0 in figure 4. They arise through the interaction of the two complexconjugate MRI harmonics with ± (ω 1 ) and the temporal flow mean. This mode is axisymmetric in ϕ (m = 0) with axial wavenumber k = 6. In the decomposition of u ϕ and B ϕ , the 0-mode corresponds to the axisymmetric mean flow, and has the largest amplitude; its spatial structure in z can be interpreted as a perturbation to the mean profile by the standing wave with ω 1 . By averaging this mode over z, we obtain the global mean profile of the DMD decomposition. Figure 6b compares this mean to the global spatio-temporal mean of the flow, and to the imposed laminar velocity profile (2.4). As MRI turbulence develops the turbulent angular momentum transport modifies the imposed rotation profile; it become flatter in the bulk and develops large gradients at the walls. This is not expected in a real astrophysical object, where turbulent fluctuations are thought to play a secondary role compared with gravity and mean rotation. The mean B φ flattens in a similar way and is also captured by DMD. has a frequency (ω 2 ) ≈ 1.1 (ω 1 ), slightly faster than the frequency of the dominant mode. It also has the same periodicity in z and ϕ, (k, m) = (3, 1), as shown in figure 5b. In the following, we will refer to the appearance of modes with comparable frequency content as mode splitting. The second mode, 2 and 2 * , has a slower temporal evolution of (ω 2 ) ∝ (ω 2 ) − (ω 1 ), and is a result of interaction between modes 1 and 2 through nonlinear terms. It has double periodicity in z, k = 6, and is axisymmetric with m = 0, as shown in figure 5c. This slow frequency, which has emerged in a secondary Hopf bifurcation, is the one forming the torus in figure 3c, and is responsible for periodic oscillations of G and flow energy E. In figure 6c, we plot the fast Fourier transform (FFT) of E mag of the flow for Ha = 145, tracking this slower oscillation in the DNS. Note that both energies and torque are quadratic quantities (2.5), so their frequencies are twice the frequencies detected in u and B. In the neighbouring figure 6d, this frequency is compared with the slow frequency detected in the DMD decomposition of B z . In the respective region of Ha, both frequencies increase when Ha decreases, with a good comparison between the two.

(c) Doubly periodic oscillation and transition to chaos
Now we consider the case of Ha = 140 (figure 4c). As discussed before, with further decrease in Ha wrinkles develop on the chaotic attractor until the flow becomes mildly chaotic (figure 3b). In this case, DMD decomposition becomes less robust but nevertheless, it is possible to detect further mode splitting in figure 4c. Now both ω 1 and ω 2 are accompanied by neighbouring modes with comparable frequencies. The interaction of these modes leads to the appearance of even slower modulations in the flow, with respective slow frequencies clustered around the mode 0, behind stronger signals from ω 2 ; this slow timescale may be the signature of the torus approaching a periodic orbit or fixed point in the flow.
In contrast to the purely periodic case with only one dominant frequency in G and E, the FFT reveals that here the signal is not perfectly doubly periodic and contains several frequency components. However, most of them can be identified as interactions between the previously detected frequency ω 2 and the modulating one with (ω 3 ) ≈ 0.01Ω i . In this regime, the former saturates at (ω 2 ) ≈ 0.045Ω i , and the modulation becomes slower as Ha decreases (figure 6d). We seek modes with a similar frequency component in the DMD decomposition of b z , which has the lowest data rank compared to the rest of the flow for all Ha, and identify a second mode splitting of (ω 3 ) ≈ 1.025 (ω 1 ) (figure 4c). The new mode has a larger axial wavelength with (k, m) = (2, 1) compared to the mode 1 (figure 5d). It was absent in more regular states of the flow and is not a result of harmonic self-interaction of the unstable modes, as its axial wavenumber is smaller. This mode is accompanied by a slow harmonic with (ω 3 ) ≈ 0.01Ω i , which is a large-scale, axisymmetric structure of k z = 1, m = 0 (figure 5d), indicating triadic interaction among the modes (ω 1 , ω 3 , ω 3 ). It is unclear whether mode 3 or 3 is of primary importance. In DMD of other flow variables, mode 3 appears more consistently than mode 3 , and tends to have a higher optimal amplitude. In figure 6d, we compare (ω 3 ) with the slow modulating frequency of E mag , and observe that the two are in agreement.

Discussion and outlook
In this work, we have employed the data-driven analysis to track transition to chaos in TCF subject to an azimuthal magnetic field. There, MRI arises as a standing wave through a supercritical Hopf bifurcation. In fluids with low conductivity (low Pm), a secondary subcritical Hopf bifurcation exists with an unstable edge state separating the periodic MRI and chaos [10]. On the contrary, a fluid with high conductivity (Pm = 1) shows more prolonged transition with diverse flow states [10]. Here, we focused on one region of this transition, Ha ∈ (120, 150), for fixed Re = 250. With decreasing Ha, the friction on the cylinders changes from constant to oscillating, and then to a modulated signal, before becoming chaotic. This transition seemingly follows a well-known Ruelle-Takens scenario, with a cascade of two Hopf bifurcations, the first at the onset of MRI, and the second when periodic oscillations of G and E develop. In the phase space, it involves a periodic orbit, a torus, and then a breakdown of the torus through folding and wrinkling of the attractor ( figure 3). On the other hand, the alignment of the attractor remains relatively unchanged, despite developing chaotic dynamics.
We employed DMD to identify the changes in the flow responsible for the temporal behaviour of its friction and energy. The first transition from the MRI standing wave (mode 1) to periodic oscillations in G happens in a process of mode splitting, i.e. mode 2 similar to the MRI mode with a slightly different frequency appears in the domain. Intuitively, the appearance of the oscillation in G and E in this case can be understood as a symmetry breaking in the system. A flow with two dynamical components, rotating in ϕ at different frequencies, is no longer invariant in ϕ in terms of the integral quantities. The friction on the cylinders at any time depends on the particular alignment of the two dynamical flow structures, periodically returning to their initial configuration.
The next flow state, where G is modulated by a slower frequency, is mildly chaotic in r-and ϕ-directions of the flow; however, it remains relatively ordered in the z-direction. The DMD decomposition of b z shows the appearance of new modes 3 and 3 with larger axial wavelengths; their footprint is visible in the rest of the flow variables. They create a frequency content comparable to the modulating frequency of the torque. The frequency (ω 3 ) is related to positive and negative regions of the field in figure 5e interchanging their location along z. As the flow becomes more complex in this case, with the modes 1, 2, 2 also influencing dynamics (figure 4c), DMD detects the slow modulation as a set of modes of similar frequency content and spatial shape, slightly different for different flow variables. Although we attribute the modulation in G, E kin and E mag to the presence of modes 3 and 3 , the modulation of G and E is likely a cumulative effect of all these harmonics. On the other side, the ω 3 modal component of magnetic field with k = 1 and k = 2 is clearly present in the dynamics of magnetic field in figure 2f and also in its space-time plots (not shown here). Thus, DMD was able to identify the flow components relevant for transition to chaotic dynamics.
As the transition to chaos proceeds further, DMD represents the flow as a set of splitting frequencies clustered about the originally dominant MRI modes (figure 4d). There are still only a few modes with large amplitudes, highlighting the low-dimensionality of this chaotic attractor. The DMD modes of chaotic flow lie inside the unit circle and appear dampened. However, their instantaneous temporal coefficients (not shown here) have chaotic rather than decaying dynamics, indicating that the linearity assumption of (3.2) is no longer valid and nonlinear dependencies between the modes emerge. This is not a concerning issue here, since we used DMD not for reduced-order modelling of the system, but as a diagnostic tool for its dynamical behaviour. The future work will include relating the temporal evolution of the modes into a nonlinear model, together with improving robustness of the presented DMD method by taking into account flow symmetries [20,21], or harnessing statistical properties of the flow [22]. Such a model could provide a quantitative description of nonlinear interactions accompanying the transition to MRI turbulence.